function [z,xc,yc] = exp_surface(b,b_range,dim,t)

br_min = min(real(b_range));
br_max = max(real(b_range));
bi_min = min(imag(b_range));
bi_max = max(imag(b_range));

br_delta = (br_max-br_min)/dim;
bi_delta = (bi_max-bi_min)/dim;

%[br,bi] = ndgrid(br_min:br_delta:br_max,bi_min:bi_delta:bi_max);
br = [br_min:br_delta:br_max];
bi = [bi_min:bi_delta:bi_max];
xc = br;
yc = bi;

z = zeros(dim+1);
s = exp(b*t);
for ii=1:dim+1,
    for jj=1:dim+1,
        ssyn = exp((br(ii)+j*bi(jj))*t);
        a = (ssyn'*s)/(ssyn'*ssyn);  % use this to find optimal amplitude as other parameters vary
        %a = 1;    % use this to fix amplitude as other parameters vary
        resid = s - a*ssyn;
        z(ii,jj) = resid'*resid;
    end
end 